###################################################
## Demographic Analysis
## SCRIPT 3: Analyze pop. tolerance tests simulated
## in script 1.
##
## For memo, creates figs/poptol_all.pdf.
## Also creates other figures we did not use.
##
## 5/17/21
###################################################

library(ppmf)
library(stringr)
library(geomander)
library(tidyverse)
library(sf)
library(redist)
library(patchwork)

la_all <- readRDS("../../data/LA/la.Rds") %>%
  rename(sld_upper = sld_up,
         sld_lower = sld_low) %>%
  filter(sld_lower != 106)

la_map <- la_all %>%
  st_set_crs(26915) %>%
  st_transform(crs = 4269) %>%
  redist_map(
    total_pop = pop,
    pop_tol = 0.005,
    existing_plan = sld_upper
  )

la_map4 <- la_all %>%
  st_set_crs(26915) %>%
  st_transform(crs = 4269) %>%
  redist_map(
    total_pop = v4_pop,
    pop_tol = 0.005,
    existing_plan = sld_upper
  )

la_map12 <- la_all %>%
  st_set_crs(26915) %>%
  st_transform(crs = 4269) %>%
  redist_map(
    total_pop = v12_pop,
    pop_tol = 0.005,
    existing_plan = sld_upper
  )

la_map19 <- la_all %>%
  st_set_crs(26915) %>%
  st_transform(crs = 4269) %>%
  redist_map(
    total_pop = v19_pop,
    pop_tol = 0.005,
    existing_plan = sld_upper
  )

## Calculate parity under existing plan
  ## creates parity_lines object to add to graph
source("xx_calculate_parity.R")

########################################################
# setup file paths
paths <- list.files("../../data/LA/sim/parity/")
types <- str_split(paths, "_") %>% do.call(rbind, .)
types <- as_tibble(types)
colnames(types) <- c("dataset", "num")
types$num <- gsub(".rds", "", types$num)

# store values here
res_df <- NULL

for (i in 1:length(paths)) {
  cat(i, "\n")
  df <- readRDS(paste("../../data/LA/sim/parity/", paths[i], sep = ""))
  df <- df %>% filter(draw != "2010 Adopted" & draw != 1)

  tb_df <- tibble(par = redist.parity(get_plans_matrix(df),
                                      la_map$pop),
                  v4_par = redist.parity(get_plans_matrix(df),
                                         la_map$v4_pop),
                  v12_par = redist.parity(get_plans_matrix(df),
                                          la_map$v12_pop),
                  v19_par = redist.parity(get_plans_matrix(df),
                                          la_map$v19_pop))

  tb_df$dataset <- types$dataset[i]
  tb_df$num     <- types$num[i]

  res_df <- bind_rows(res_df, tb_df)

}

# for theme_ppmf()
source("../00_custom_functions.R")

## Creates figs/poptol_all.pdf for main text
p_line <- res_df %>%
  filter(num != 0.5 & dataset != "smc19") %>%
  mutate(num = as.numeric(num)) %>%
  pivot_longer(cols = par:v12_par) %>%
  group_by(dataset, num, name) %>%
  summarise(invalid_pct = mean(value > num)) %>%
  mutate(name = case_when(name == "par" ~ "Census",
                          name == "v4_par" ~ "DAS-4.5",
                          name == "v12_par" ~ "DAS-12.2"),
         dataset = case_when(dataset == "smc" ~ "Sampled from: Census",
                             dataset == "smc4" ~ "Sampled from: DAS-4.5",
                             dataset == "smc12" ~ "Sampled from: DAS-12.2")) %>%
  arrange(dataset, desc(num)) %>%
  ggplot(aes(x = num,
             y = invalid_pct,
             color = name, group = name)) +
  geom_line() + geom_point(size = 1.5) +
  geom_vline(data = parity_lines %>% filter(name != "DAS-19.61"),
             aes(xintercept = int),
             lty = "dashed", col = "black", lwd = 0.2) +
  geom_text(data = parity_lines %>% filter(name != "DAS-19.61"), inherit.aes = FALSE,
            aes(label = paste("Enacted: \n",
                              round(int * 100, 2), "%", sep = "")),
            x = -0.85, y = 0.8,
            size = 2.5, family = "Times",
            color = c("#29483A", "#D94154", "#D94154")) +
  facet_wrap(~dataset) +
  ## color scheme
  ## from https://github.com/CoryMcCartan/wacolors/blob/main/R/colors.R
  scale_color_manual(name = "Evaluation\ndata source:",
                     values = c("#9FB6DA", "#759C44", "#29483A")) +
  labs(x = "Intended Population Tolerance",
       y = "Fraction of Plans Invalid") +
  scale_x_log10(breaks = c(0.001, 0.005, 0.01, 0.05,
                           0.1, 0.2, 0.5),
                labels = c("0.1%", "0.5%", "1%", "5%",
                           "", "20%", "50%")) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  theme(strip.text.x = element_text(size = 8, face = "bold"),
        legend.title = element_text(size = 9),
        legend.text = element_text(size = 8),
        axis.text = element_text(size = 5.5),
        axis.title = element_text(size = 8)) +
  theme_ppmf()
## ggsave("../../figs/poptol_all.pdf", p_line, height = 2, width = 6.5)

## postscript plot
p_line_post <- res_df %>%
  filter(num != 0.5 & dataset %in% c("smc", "smc19")) %>%
  mutate(num = as.numeric(num)) %>%
  pivot_longer(cols = c(par, v19_par)) %>%
  select(-v4_par, -v12_par) %>% 
  group_by(dataset, num, name) %>%
  summarise(invalid_pct = mean(value > num)) %>%
  mutate(name = case_when(name == "par" ~ "Census",
                          name == "v19_par" ~ "DAS-19.61"),
         dataset = case_when(dataset == "smc" ~ "Sampled from: Census",
                             dataset == "smc19" ~ "Sampled from: DAS-19.61")) %>%
  arrange(dataset, desc(num)) %>%
  ggplot(aes(x = num,
             y = invalid_pct,
             color = name, group = name)) +
  geom_line() + geom_point(size = 1.5) +
  geom_vline(data = parity_lines %>% filter(name %in% c("Census", "DAS-19.61")),
             aes(xintercept = int),
             lty = "dashed", col = "black", lwd = 0.2) +
  geom_text(data = parity_lines %>% filter(name %in% c("Census", "DAS-19.61")), 
            inherit.aes = FALSE,
            aes(label = paste("Enacted: \n",
                              round(int * 100, 2), "%", sep = "")),
            x = -0.85, y = 0.8,
            size = 2.5, family = "Times",
            color = c("#29483A", "#29483A")) +
  facet_wrap(~dataset) +
  ## color scheme
  ## from https://github.com/CoryMcCartan/wacolors/blob/main/R/colors.R
  scale_color_manual(name = "Evaluation\ndata source:",
                     values = c("#9FB6DA", "#465177")) +
  labs(x = "Intended Population Tolerance",
       y = "Fraction of Plans Invalid") +
  scale_x_log10(breaks = c(0.001, 0.005, 0.01, 0.05,
                           0.1, 0.2, 0.5),
                labels = c("0.1%", "0.5%", "1%", "5%",
                           "", "20%", "50%")) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  theme(strip.text.x = element_text(size = 8, face = "bold"),
        legend.title = element_text(size = 9),
        legend.text = element_text(size = 8),
        axis.text = element_text(size = 5.5),
        axis.title = element_text(size = 8)) +
  theme_ppmf()
## ggsave("../../figs/poptol_all_post.pdf", p_line_post, height = 2, width = 6.5)
